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ABSTRACT 


A  novel  technique  to  determine  the  phase  velocity  of  long-wavelength 
shoaling  waves  is  investigated.  Operationally,  the  technique  consists  of  three 
steps.  First,  using  the  Hilbert  transform  of  a  time  series,  the  phase  of  the 
analytic  signal  is  determined.  Second,  the  correlations  of  the  phases  of  analytic 
signals  between  two  points  in  space  are  calculated  and  an  average  time  of  travel 
of  the  wave  fronts  is  obtained.  Third,  if  directional  spectra  are  available  or  can 
be  determined  from  time  series  of  large  array  of  buoys,  the  angular  information 
can  be  used  to  determine  the  true  time  of  travel.  The  phase  velocity  is  obtained 
by  dividing  the  distance  between  buoys  by  the  correlation  time.  Using  the  Hilbert 
transform  approach,  there  is  no  explicit  assumption  of  the  relation  between 
frequency  and  wavenumber  of  waves  in  the  wave  field,  indicating  that  it  may  be 
applicable  to  arbitrary  wave  fields,  both  linear  and  nonlinear.  Limitations  of  the 
approach  are  discussed. 
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I.  INTRODUCTION 


As  ocean  waves  shoal,  the  wave  field  evolves  substantially  from  its  deep 
water  state.  While  the  geometry  of  the  dynamic  currents  and  of  the  bottom 
topography  determine  the  direction  and  to  some  extent  the  energy  of  the  spectral 
components,  nonlinearities  in  the  wave  field  lead  to  harmonic  generation  and  to 
spectrum  broadening.  Furthermore,  for  a  nonlinear  wave  field,  there  is  an 
apparent  negative  refraction  of  the  spectral  peak  (Abreu  et  al.  1992).  Thus,  for  a 
known  sloping  bottom,  the  refraction  away  from  the  beach  normal  is  a  direct 
consequence  of  nonlinear  energy  transfer.  On  the  other  hand,  linear  theory 
does  not  predict  most  of  these  changes,  and  for  a  given  shoaling  spectrum, 
linear  theory  might  lead  to  an  incorrect  determination  of  the  beach  normal. 

Directional  spectrum  has  become  commonly  used  to  describe  nonlinear 
random  ocean  waves  because  it  provides  ways  to  estimate  the  significant  wave 
heights  in  a  breaker  and  also  to  determine  the  predominant  direction  for  net 
sediment  transport.  Considerable  progress  has  been  made  over  the  last  two 
decades  to  study  the  evolution  of  the  wave  spectrum  due  to  refraction, 
diffraction,  and  nonlinear  interaction  under  conditions  of  mild  slope  (Lui  et  al. 
1985,  Freilich  and  Guza  1984,  Abreu  et  al.  1992).  However,  field  conditions 
often  violate  the  basic  assumption  that  the  wavelength  be  smaller  than  the 
characteristic  length  scale  of  bottom  variations.  Nevertheless,  it  can  be  shown 
(Larraza,  1994)  that  using  a  ray  optics  approximation,  the  determination  of  the 
kinematic  wave  parameters  (wavelength,  direction  of  propagation)  for  a  given 
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geometry  is  within  5%  of  its  true  value.  On  the  other  hand,  the  determination  of 
the  dynamic  values  (energy,  momentum)  exceeds  at  least  15%  the  true  value. 
Thus  by  using  ray  theory,  sediment  transport  and  other  relevant  effects  caused 
by  the  wave's  momentum  will  be  ill  estimated. 

For  an  arbitrary  bottom  topography,  wave  energy  and  momentum  can  be 
calculated  within  the  physical  optics  approximation.  The  passage  from 
geometrical  to  physical  optics  can  in  principle  be  accomplished  by  use  of  a  path 
integral  formulation  approach.  Within  this  formulation,  the  basis  states  are 
characterized  by  the  "rays",  and  the  physical  optics  results  from  a  sum  over  all 
the  possible  paths  (and  not  only  the  path  determined  by  Fermat's  principle).  The 
path  integral  formulation  approach  has  not  been  considered  previously  either  for 
gravity  waves  over  a  changing  topography  or  for  shoaling  waves.  This  situation 
is  particularly  surprising  because  the  problem  can  be  stated  from  the  onset  by 
using  the  known  results  of  ray  optics,  either  linear  or  nonlinear,  and  performing  a 
sum  over  all  possible  ray  paths  of  the  basis  states.  Even  though  there  are 
standard  analytical  techniques  to  deal  with  certain  special  cases,  the  most 
general  ones  can  only  be  solved  numerically. 

Because  the  sea  floor  is  known  to  be  a  controlling  factor  in  low-frequency 
shallow  water  wave  propagation  (Akal,  1980,  Akal  and  Jensen  1983),  phase 
velocity  measurements  can  in  principle  be  used  to  determine  bottom  topography, 
which  can  in  turn  be  used  for  true  ground  determination  in  remote  sensing 
applications.  The  inverse  problem,  that  is,  the  determination  of  the  topography 
given  the  evolution  of  the  spectrum  along  a  path  requires  an  accurate  technique 
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to  determine  celerity  and  direction  of  the  waves  in  a  given  wave  field  data. 
Common  techniques  of  Fourier  transform  are  usually  inadequate.  In  the  Fourier 
transform,  a  real  space-time  signal  is  converted  to  a  complex  frequency- 
wavevector  signal.  Thus,  for  applications  to  ocean  wave  field  data  a  dispersion 
relation  between  frequency  and  wavenumber  has  to  be  assumed  in  order  to 
determine  the  phase  velocity  of  the  spectral  components.  The  problem  can  be 
complicated  by  nonlinearities  in  the  wave  field. 

On  the  other  hand,  the  Hilbert  transform  of  a  real  valued  time-space  signal 
is  another  real  value  time-space  signal  for  which  the  phase  of  a  signal  can  be 
calculated  in  terms  of  the  signal  itself  and  Its  Hilbert  transform.  Thus  by  using 
the  time  series  of  a  signal  measured  at  different  points,  one  can  determine  a 
relative  phase  shift  and  extract  phase  velocity  measurements  without  reference 
to  a  dispersion  relation. 

In  this  thesis  we  present  an  analysis  of  ocean  surface  waves  data 
obtained  by  The  Naval  Research  Laboratory  Stennis  Center,  in  support  of 
Hamlet’s  Cove  I  (Smith,  1994).  An  immediate  goal  is  to  extract  phase  velocity 
information  from  time  series  at  three  buoys.  The  data  analysis  technique  used 
throughout  this  thesis  is  based  on  the  Hilbert  transform.  This  technique  allow  us 
to  get  local  properties,  local  energy  and  local  frequency,  of  the  signal  instead  of 
the  global  properties  that  others  techniques,  like  the  Fourier  transform,  would 
give  (Long,  1995). 

In  Chapter  II  we  give  a  brief  introduction  to  linear  water  waves  and  point 
out  their  dispersive  character.  We  also  present  a  description  of  the  local 
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topography  for  the  wave  data  that  we  analyze  and  emphasize  the  need  for  an 
inverse  problem  in  arbitrary  topographies.  Chapter  III  gives  a  detailed  account 
of  the  basic  ideas  involved  in  Hilbert  transforms.  In  Chapter  IV  we  present  the 
results  and  the  procedures  used  for  the  analysis  of  the  data.  The  data  was 
analyzed  using  the  Fast  Fourier  Transform  (FFT),  the  Hilbert  transform,  and  the 
cross-correlation  functions  of  MATLAB.  In  Chapter  V  we  discuss  the  results  and 
provide  some  recommendations  for  future  work. 
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II.  SURFACE  WAVES  IN  THE  HAMLET’S  COVE  I  EXPERIMENT 


In  the  first  part  of  this  chapter,  we  present  a  brief  introduction  to  surface 
gravity  waves  over  uniform  depth  and  emphasize  their  dispersive  behavior.  In 
contrast,  the  data  that  we  analyze  throughout  this  thesis  is  over  nonuniform 
terrain.  The  second  part  of  this  chapter  gives  a  detailed  description  of  the  local 
topography  for  the  Hamlet’s  Cove  I  experiment. 

A.  SURFACE  WAVES  IN  WATER 

When  a  group  of  waves  moves  across  the  surface  of  water,  each 
particular  wavecrest  travels  faster  than  the  group  as  a  whole,  and  eventually 
passes  through  it.  Thus  new  crest  continually  are  being  created  at  the  back  of 
the  group  while  old  crests  are  disappearing  at  the  front. 

This  behavior  does  occur  because  water  waves  are  dispersive.  We  can 
say  that  this  implies  that  the  different  Fourier  components  that  make  the  general 
disturbance  have  phase  velocities  that  depend  on  their  wavelength.  For  surface 
waves,  if  the  surface  elevation  of  water  with  uniform  depth  h  is  described  by  the 
wavetrain 

ri  =  A  cos  (kx-  cot)  ,  (1.1) 

the  wave  speed  is  given  by 
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c  =  CO  /  k  = 


7-  tanh(kh) 

U 


n1/2 

/ 


(1.2) 


so  that  waves  of  longer  wavelength,  A,  =  2k  /  k,  travel  faster.  In  Eq  (1 .1)  k  is  the 
wavenumber  and  co  is  the  frequency  of  the  wave. 

While  each  individual  wavecrest  travels  with  speed  c,  the  velocity  of  travel 
of  the  group  as  a  whole  is  the  group  velocity 

Cg  =  dco/dk  (1.3) 

The  phase  velocity  in  water  with  uniform  depth  h  (Acheson,  1995)  can  not 
exceed  ■  For  kh  large,  i.e.  h  »X,  tanh(k/7)  =  1  ,  and 

c^  =  g  /  k  ,  (1.4) 

corresponding  to  the  case  of  infinite  depth.  A  good  approximation  for  deep  water 
waves  is  for  depths  h  greater  than  about  In  shallow  water,  i.e.  h«XI  2k, 

tanh(kh)  =  kh  ,  and  equation  (1 .2)  becomes 

c^  s  gh  .  (1.5) 
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This  equation  is  important,  because  as  we  can  see  c  is  independent  of  k  for 
shallow  water.  As  we  will  see,  the  data  analyzed  in  this  thesis  is  remarkably 
nondispersive. 

B.  HAMLET’S  COVE  TOPOGRAPHY 

The  data  analyzed  in  this  thesis  is  from  three  bottom-mounted  pressure 
gauge/  current  meters  for  the  measurements  correspond  to  directional  wave 
spectra,  currents,  and  tidal  elevation.  The  gauges  were  deployed  in  depths  of 
10’ ,  15’,  and  25’  installed  in  support  of  Hamlet  Cove  I  experiment. 

The  three  bottom-mounted  pressure  gauges  were  placed  on  either  side  of 
the  sand  bar  as  shown  in  Figures  2.1  and  2.2 

86M8’10”  86M8’20” 

30°23’10” 

30«23’00” 

30°22’50” 


•  #3  (Depth  10’) 

#2  (Depth  15’) 

#  1  (Depth  25’r 

Figure  2.1 .  Coordinates  for  the  Positions  of  the  gauges  and  their  depths  in  feet. 
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Figure  2.2.  Distance  from  shore  (m)  of  the  three  buoys  and  reference  to  the  sand  bar.  The  sand 
bar  profile,  measured  relative  to  the  surface,  is  a  t  a  depth  of  1 .3  m  1 50  m  off  shore. 

The  buoys  used  collected  data  from  the  wave  as  pressure,  x-component 
of  the  velocity  and  the  y-component  as  well.  The  sampling  was  continuous  with 
a  sampling  frequency  of  2  Hz  . 

Because  the  pressure  gauge  responds  to  the  rise  and  fall  to  the  free 
surface,  the  pressure  at  depth  z  attenuates  by  a  factor  proportional  to  e  ' 

where  k  is  2TdX  and  %  is  the  wavelength  of  the  Fourier  component  of  the  surface 
disturbance.  At  the  15’  and  25’  depths,  (gauges  1  and  2)  frequencies  above  0.3 
Hz  were  below  the  noise  level  of  the  instruments  and,  therefore,  were  ignored  in 
the  analysis.  The  high  frequency  cut-off  for  gauge  number  3  is  0.5  Hz. 

Table  2-1  shows  the  data  collection  periods  for  each  instrument.  We  analyzed 
the  data  collected  by  the  three  buoys  in  16  July  1994. 
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Gauge 

1 

2 

3 

DATES 

14JUL-20  JUL 

14  JUL-  20  JUL 

14  JUL-  20  JUL 

20  JUL-29  JUL 

20  JUL-  29  JUL 

20  JUL-  29  JUL 

29  JUL-1  AUG 

29  JUL-  3  AUG 

29  JUL-  3  AUG 

3  AUG-10  AUG 

3  AUG-10  AUG 

3  AUG-9  AUG 

12  AUG-17  AUG 

10  AUG-17  AUG 

FAILED 

18  AUG-24  AUG 

18  AUG-24  AUG 

17AUG-24  AUG 

24  AUG-4  SEP 

FAILED 

FAILED 

Table  2-1  Dates  for  the  data  collected  at  the  three  buoys.  The  data  analyzed  in  this  thesis 
corresponds  to  data  collected  on  16  July  1994. 


From  the  bathymetry  survey  shown  in  Figure  2.2,  we  can  see  that 
refraction  and  diffraction  effects  from  waves  60  m  long  and  longer  may  be  used 
to  characterize  the  bottom  topography.  This  long  waves,  due  mainly  to  swell,  are 
essentially  nondispersive  and  already  violate  the  conditions  for  ray  optics  to 
apply.  Because  shorter  dispersive  waves  are  mainly  due  to  local  winds,  they 
would  not  give  accurate  topography  estimates.  Thus,  phase  velocity  estimates 
from  the  long  waves  are  the  first  step  to  determine  an  unknown  bathymetry. 
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III.  THE  HILBERT  TRANSFORM 

In  this  chapter  we  review  some  of  the  basic  concepts  of  the  Hilbert 
transform,  following  the  review  article  by  Bendat  (1989).  In  order  to  elucidate  the 
concept  and  possible  applications  of  the  Hilbert  transform,  we  present  three 
equivalent  definitions,  namely  as  a  convolution  integral,  as  a  7i/2  phase  shifter, 
and  as  an  imaginary  part  of  an  analytic  signal. 

The  Hilbert  Transform  of  a  real  valued  signal  is  a  complex  signal  called 
the  analytic  signal.  The  real  part  of  the  analytic  signal  is  the  original  real-valued 
time  signal  while  the  imaginary  part  is  a  copy  of  the  original  signal  with  each  of  its 
Fourier  components  shifted  in  phase  by  90  (Bendat,  1989) .  The  transformed 
signal  retains  the  same  amplitude  and  frequency  information  contained  in  the 
original  signal  with  the  added  phase  information  dependent  on  the  phase  of  the 
original  data. 

For  any  real  time  series,  ^(t),  the  Hilbert  transform  ^  (t)  is 


in  which  P  implies  the  principal  value.  The  analytic  continuation  of  the  Hilbert 
transform  is  thus  the  Cauchy  integral  which  corresponds  to  translations  of  an 
analytic  function  in  the  complex  domain  of  integration. 

The  analytic  signal  Z(t)  is  the  real  based  time  series  defined  by 
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z  (t)  =  C(t)  +  j  ^(t) 


(3.2) 


which  can  also  be  written  as 

Z(t)  =  A(t)  e  ,  (3.3) 

where  A(t)  is  called  the  envelope  signal  or  the  amplitude  and  i}(t)  is  called  the 
instantaneous  phase  signal,  defined  by: 


A(t)  =  [C"(t)  + 


and 


m  =  tan-^[4(t)  /  C(t)]. 


respectively.  The  instantaneous  frequency  co  is  given  by 


CO  = 


d^  (t) 
dt 


(3.4) 


(3.5) 


(3.6) 


and  the  local  energy  is  one  half  of  the  square  value  of  the  amplitude.  Eqs.  (3.4), 
(3.5),  and  (3.6)  represents  local  properties  in  time. 

The  properties  of  the  Hilbert  Transform  can  be  best  understood  in  terms 
of  three  equivalent  definitions. 

1)  Definition  as  Convolution  Integral 

The  Hilbert  Transform  of  a  real  valued  function  ^(t)  extending  from  -  <==  <  t 
<  +  «=  is  a  real  valued  function  as  described  in  (3.1): 


12 


(3.7) 


thus  ^  (t)  is  the  convolution  integral  of  ^(t)  with  ( 1  /  rt  t),  written  as 

Ut)  =  ?  (t)  *(-^)  .  (3.8) 

V7l  Ij 

2)  Definition  as  (tt  /  2)  Phase  Shift  System 

If  F^(  f )  is  the  Fourier  transform  of  ^  (t) ,  namely 

+00 

F  ^(f)  =  (t)  e  dt  ,  (3.9) 

—DO 


then,  from  the  convolution  property  of  the  Hilbert  transform  (Bendat,  1989)  and 
equation  (3.7),  it  follows  that  F^(  f )  is  the  Fourier  transform  of  ^(t),  multiplied  by 
the  Fourier  transform  of  (1  /  ti  t)  where 


J 


e 

- ; — dt  =  -jsgnf  = 

TCt 


- j  for  f  >  0 
j  for  f  <  0 


(3.10) 


Hence,  equation  (3.6)  is  equivalent  to  the  passage  of  ^  (t)  through  a  system 
defined  by  ( -j  sgnf)  to  yield: 

F'(f)  =  (-jsgnf)F''(f)  ,  (3.11) 

where  F^^(f)  is  the  Fourier  transform  of  ^(t).  The  complex-valued  quantity  F^(  f )  is 
the  Hilbert  transform  of  the  complex-valued  quantity  F^^  ( f )  as  defined  by: 
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F'(f)  =  H[F"(f)]=  (-isgnf)F"(f) 


(3.12) 


The  Fourier  transform  ( -j  sgnf)  can  be  represented  by 


B  (f)  =  -  j  sgnf  = 


e  ’  ^  ^Mor  f  >  0 
e  <  0 


(3.13) 


or  using  the  complex  polar  notation 


B(f)  =  |B(f)|  e  “  . 


(3.14) 


Hence,  B(f)  is  a  [n!  2)  phase  shift  system  where 


^(f)|=  1  for  all  f , 


(3.15) 


and 


<Pb(f) 


j  7c  /  2  for  f  >  0  1 
[-  7t  /  2  for  f  <  0  I  ’ 


(3.16) 


Thus  if  one  writes: 


(3.17) 


it  foiiows  that 


F'(f)  =  |F'(f)|  =  |F'^(f)|  e''l ^ 


(3.18) 
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Thus  the  Hilbert  transform  consist  of  passing  ^(t)  through  a  system  which  leaves 


the  magnitude  of  v"  (f)  unchanged,  but  changes  the  phase  from  ([)  x  (f)  to  ([)  x  (f)  + 
(j)  b  (f).  Thus 

<Px(f)  9x  W  +  (tt  /  2)  for  f  >  0  ,  (3.19a) 

9x(f)  ^  9x  (f)  ■  (^  /  2)  for  f  <  0  (3.19b) 

in  other  words,  the  Hilbert  transform  shifts  by  tc  /  2  for  positive  frequencies  and 
by  -Td2  for  negative  frequencies. 

3)  Definition  as  imaginary  part  of  Anaiytic  Signai 

A  third  useful  way  to  understand  and  to  compute  the  Hilbert  transform  is 
to  introduce  the  analytic  signal  Z(t)  which  allows  us  to  work  with  the  amplitude 
A(t),  and  the  phase  <])  (t)  equations  (3.4  and  3.5).  Time  derivatives  of  the  phase 
yield  the  instantaneous  frequency  co  (3.6). 

Examples  of  Hilbert  transforms  and  the  corresponding  Fourier  transform 
of  some  simple  functions  are  shown  in  Figure  3.1 .  The  convolution  and  phase 
shifts  due  to  the  Hilbert  transform  are  apparent.  Also  apparent  is  the  fact  that 
the  analytic  signal  preserves  local  amplitude. 


SIGNAL  TRANSFORM 

HILBERT  TRANSFORM 

FOURIER 
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IV.  DATA  ANALYSIS 


In  this  chapter  we  present  the  procedure  and  the  results  of  the  analysis  of 
the  Hamlet’s  Cove  I  surface  wave  data.  The  main  goal  is  to  determine  if  there  is 
any  phase  correlation  between  data  at  the  three  buoys  after  low  pass  filtering 
and  applying  the  Hilbert  transform  to  the  filtered  data.  Phase  correlations  can 
provide  valuable  information  because  they  can  tell  us  if  the  signal  recorded  from 
the  three  buoys  corresponded  to  the  same  wave  and  if  so,  they  can  yield  directly 
the  phase  velocity  of  the  wave.  Phase  velocity  is  obtained  by  interpolating  the 
correlation  time  between  the  buoys  for  which  the  location  is  known. 

The  data  analyzed  is  of  pressure  measurements  from  a  pressure  gauge 
and  two  orthogonal  components  of  velocity  from  a  current  meter.  Because 
pressure  gauge  measurements  are  more  sensitive  to  long  waves,  we  will 
consider  the  low  pass  filtered  output  from  for  pressure  gauge  of  measurements 
on  16  July  1994. 

A.  THE  CHEBYSHEV  TYPE  II  FILTER 

An  optimal  filter  needs  to  maintain  amplitude  information  with  no  phase 
distortion  in  the  pass  band.  In  our  case,  the  latter  requirement  is  particularly 
important  because  from  our  approach,  phase  velocity  values  can  only  be 
determined  from  phase  correlations.  Because  the  filtering  process  should 
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minimize  distortion  of  all  small  amplitude  fluctuations,  the  stop  band  has  to  be  at 
least  40  dB  below  the  pass  band  to  eliminate  possible  leakage  problems  and  it 
has  to  have  a  minimum  transition  width  to  ensure  good  frequency  separation. 
Most  of  these  requirements  can  be  met  by  with  a  Chebyshev  type  II  filter.  Higher 
order  filters  could  introduce  numerical  errors,  so  the  filter  has  to  meet  the 
requirements  with  the  lowest  order  possible. 

As  shown  in  the  Appendix  (mfile),  we  used  Matlab  Signal  Processing 
Toolbox  (Mathworks,  1994).  Of  the  two  techniques.  Analog  Prototyping  and 
Direct  Design,  we  chose  Analog  Prototyping  because  the  Direct  Design 
technique  has  very  stringent  numerical  accuracy  constraints. 

The  Analog  Prototyping  technique  allows  the  construction  of  digital 
equivalents  of  certain  classical  analog  filters,  especially  Butterworth,  Bessel, 
Chebyshev,  and  Elliptic.  Among  them  and  for  ours  purpose  the  Chebyshev  Type 
1 1  filter  has  the  narrower  transition  and  the  flatter  pass  band  response  needed  for 
the  analysis  of  the  data. 

Analog  filters  are  HR  (infinite-duration  impulse)  type  filters,  which 
necessarily  introduce  phase  distortion.  In  order  to  overcome  this  problem  the 
data  was  filtered  both  directions,  forward  and  reverse  by  using  the  Matlab 
function  filtfilt  (Mathworks,  1994).  With  this  approach,  we  can  get  precise  zero- 
phase  distortion  and  doubling  of  the  filter  order.  Also,  filtfilt  minimizes  startup 
and  ending  transients  by  adjusting  initial  conditions. 


18 


B.  THE  XCORRELATION  FUNCTION 


Matlab’s  function  xcorr  estimates  the  cross-correlation  sequence  of 
random  processes.  For  the  purpose  of  determining  phase  velocity  values  from 
apparently  random  time  series,  the  cross-correlation  function  is  important 
because  it  can  provide  a  measure  of  the  similarity  between  two  signals. 

The  xcorrfunction  can  calculate  the  cross-correlation  between  two  arrays 
of  vectors  of  the  same  length,  the  autocorrelation  for  a  vector,  and  the  cross¬ 
correlation  for  a  matrix.  In  all  these  forms  of  xcorr,  the  zero  lag  of  the  output  is  in 
the  middle  of  the  sequence.  In  order  to  normalized  the  sequence  we  used  the 
option  'coeff  so  the  autocorrelations  at  zero  lag  are  identically  1 .0. 

C.  PROCEDURE  AND  RESULTS 

Figure  4.1  shows  the  Fourier  transform  of  the  pressure  time  series  for  the 
three  buoys.  From  the  shape  of  the  spectrum,  we  can  conclude  that  waves  due 
to  the  local  wind  may  correspond  to  frequencies  above  0.1  Hz.  If  phase  velocity 
values  can  be  used  to  determine  local  bathymetry,  waves  generated  by  the  local 
wind  are  unwelcome  noise  for  the  purpose  of  our  analysis.  Instead,  the  longer 
waves  in  the  swell  may  characterize  the  bathymetry  more  accurately  because  of 
their  direct  response  to  the  bottom  topography  manifested  by  refraction  and 
diffraction  effects.  Swell  may  correspond  to  frequencies  below  0.08  Hz,  which 
we  chose  as  the  cutoff  frequency  for  the  low  pass  band  of  the  filter. 
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Figure  4.1.  FFT  of  the  pressure  data  for  the  three  buoys.  Frequencies  above  0.1  Hz  may 

correspond  to  waves  generated  by  the  local  wind  as  evidenced  by  the  increase  of 
spectral  energy  and  range. 

With  this  information  in  mind,  we  filtered  the  whole  time  series  using  a 
Chebyshev  Type  II  filter  with  a  cut  off  frequency  of  0.08  Hz.  Because  of  the 
Chebyshev  Type  II  filter  introduces  phase  distortion,  the  wave  data  was  filtered 
in  both  the  forward  and  reverse  directions  using  the  MATLAB  function  filtfilt,  thus 
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obtaining  precise  zero-phase  distortion  and  doubling  of  the  filter  order.  As 
mentioned  before,  ////fZ/f  minimizes  startup  and  ending  transients  by  adjusting 
initial  conditions. 

We  applied  the  Hilbert  transform  to  the  filtered  time  series.  Figure  4.2 
shows  the  Matiab  Power  Spectrum  function  (PSD)  of  unfiltered  time  series  and 
the  PSD  of  the  Hilbert  transform  of  the  data  from  the  first  buoy  after  fiitering. 
Apparent  from  the  figure,  is  that  the  power  spectrum  of  the  unfiltered  data  does 
not  exhibit  very  good  frequency-energy  decomposition.  On  the  other  hand,  the 
power  spectrum  of  the  Hilbert  transform  of  the  filtered  time  series  shows  good 
frequency-energy  decomposition  and  enhancement  of  the  more  energetic  waves 
of  the  fiitered  time  series. 

The  data  from  each  buoy  has  172,800  points  corresponding  to  two  data 
point  per  second  in  a  full  day.  Because  this  large  array  is  numerically  very 
demanding,  we  chose  to  work  with  a  window  of  two  hours,  or  10%  of  the  whole 
data.  We  applied  Matlab’s  Hilbert  transform  (^(t))  function  to  the  filtered  data  in 
this  time  window.  Along  with  the  real  time  series  (Y1),  MATLAB  generates  the 
analytical  signal 

Z(t)  =  Y1(t)  +  i^  (t)  ,  (4.1) 

which  is  similar  to  Eq.  (3.2). 
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Figure  4.2.  aj  Power  spectrum  density  for  the  pressure  data  at  buoy  1  without  filtering,  b)  Power 
spectra  density  of  the  Hilbert  transform  after  low  pass  filtering  using  a  Chebyshev 
type  II  filter  with  a  0.08  cutoff  frequency. 


Matlab’s  functions  angle  and  unwrap  (phase  in  radians  2*pi  wrap),  provide 
the  phase  of  the  time  series  (4.1).  Shown  in  Figure  4.3  are  the  graphs 
corresponding  to  the  phase  of  the  whole  pressure  data  from  buoy  1  and  also  for 
the  two  hours  window. 
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Figure  4.3.  Phase  (rad)  as  a  function  of  time  (min)  for  the  pressure  data  at  buoy  1  for  a)  one  day, 
b)  two  hours. 

The  local  frequency  is  by  definition  the  slope  of  the  phase  and  we  applied 
to  the  phase  data  the  MATLAB  function  diff  and  divided  by  the  time  between 
points  (called  tbp  in  Appendix)  to  get  this  slope.  Figure  4.4  shows  a  plot  of  the 
local  frequency  as  a  function  of  time  obtained  from  the  data  shown  in  Figure  4.3. 
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Figure  4.4.  Local  angular  frequency  as  a  function  of  time  for  buoy  1  for  a  two  hour  time  window. 
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Figure  4.5  shows  the  cross-correlation  of  data  from  buoys  1 , 2,  and  3 
using  the  xcorr function  in  MATLAB.  The  cross-correlation  is  blown-up  in  Figure 
4.6.  At  the  maximum  of  the  cross-correlation  we  can  obtain  the  time  when  the 
signals  from  two  different  buoys  are  highly  correlated.  Using  this  time,  and 
knowing  the  distance  between  buoys,  we  can  estimate  the  phase  velocity  at 
which  the  wave  is  moving  from  one  buoy  to  the  other.  Assuming  that  the  three 
buoys  started  to  record  the  data  at  the  same  time,  we  get  the  phase  velocities 
estimates  shown  in  Table  4-1 . 


Figure  4.5.  Cross  correlation  of  the  phase  between  a)  buoys  1  and  2,  b)  buoys  1  and  3,  and  c) 
buoys  3  and  2. 
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Table  4-1.  Relations  of  the  distances  among  the  buoys,  correlation  time  and  phase  velocity. 


25 


26 


V.  DISCUSSION  OF  RESULTS  AND  RECOMMENDATIONS  FOR 

FUTURE  WORK 


We  have  investigated  the  Hilbert  transform  as  a  tool  to  determine  phase 
velocity  from  ocean  wave  data  of  three  buoys.  The  values  tabulated  in  Table  4-1 
were  obtained  by  dividing  the  distance  between  buoys  by  the  time  for  maximum 
phase  correlation  of  the  time  series  of  pressure  measurements  at  the 
corresponding  buoys. 

Table  5-1  is  a  comparison  between  the  values  of  phase  velocity  as 
determined  by  the  correlation  times  to  the  values  of  phase  velocity  assuming 
linear  shallow  water  theory  over  uniform  bottom  topography.  The  value  of  the 
phase  velocity  determined  from  shallow  water  theory  represents  an  upper  bound. 
If  the  normal  to  the  wave  front  makes  an  angle  with  the  line  joining  two  buoys, 
the  value  of  phase  velocity  obtained  from  correlations  of  the  time  series  will  yield 
an  underestimation.  Thus,  the  low  values  obtained  from  the  correlations  may  be 
due,  in  part,  to  an  effective  average  over  all  possible  directions. 


phase  velocity  from 
correlations 

average  depth 

calculated  phase 
velocity 

Buoy  1  -  Buoy  2 

3.63  m/s 

6.10  m 

7.73  m/sec 

Buoy  2  -  Buoy  3 

1 .35  m/s 

3.81  m 

6.11  m/sec 

Buoy  1  -  Buoy  3 

3.36  m/s 

5.34  m 

7.23  m  /  sec 

Table  5.1  Comparison  of  phase  velocity  values  determined  from  correlations  (second 

column)  and  from  shallow  water  theory  (fourth  column)  assuming  uniform  depth 
(third  column). 
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It  is  then  apparent  that  to  accurately  determine  the  phase  velocity  from 
time  series,  we  require  an  accurate  directional  spectral  data  or  surface  height 
time  series  from  a  bigger  array  of  buoys.  An  accurate  directional  spectra 
provides  the  necessary  angular  information  that  can  be  used  to  determine  the 
true  value  of  the  phase  velocity.  As  shown  by  Rikiishi  (1978)  time  series  from  a 
minimum  of  three  buoys  cannot  be  used  to  determine  directional  spectra  for  a 
broadband  distribution  of  waves. 

Any  data  analysis  on  time  series  from  pressure  gauge  measurements 
have  fundamental  limitations.  First,  because  the  pressure  gauge  responds  to 
the  rise  and  fall  to  the  free  surface,  the  pressure  at  depth  z  attenuates  by  a  factor 
proportional  to  e’'^^  where  k  is  ItilX  and  X  is  the  wavelength  of  the  Fourier 
component  of  the  surface  disturbance.  Unless  the  wave  data  corresponds 
exclusively  to  shallow  water,  the  correlation  time  determined  from  the  phase  of 
the  Hilbert  transform  at  two  points  may  be  inaccurate.  Second,  if  the  wave  field 
is  nonlinear,  the  pressure  and  the  surface  height  are  nonlinearly  related.  Due  to 
nonlinearities,  the  correlation  time  from  pressure  gauge  measurements  may  yield 
misleading  phase  velocity  results.  Thus,  phase  velocity  values  from  phase 
correlations  using  the  Hilbert  transform,  are  more  reliable  with  time  series  of 
surface  height  measurements. 

If  the  deficiencies  noted  above  can  be  corrected,  the  preliminary  results  of 
this  thesis  are  evidence  of  the  powerful  technique  of  using  Hilbert  transform  to 
determine  the  phase  velocity  of  low  frequency  shoaling  waves.  This  technique 
can  be  summarize  operationally  as  follows.  The  analytic  signal  of  a  time  series 
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possesses  a  natural  phase.  Correlations  of  the  phase  yield  an  average  time  of 
travel.  If  directional  spectra  are  available,  the  angular  information  can  then  yield 
the  true  time  of  travel  from  which  the  phase  velocity  can  be  determined.  In  this 
approach,  there  is  no  explicit  assumption  of  the  relation  between  frequency  and 
wavenumber  of  waves  in  the  wave  field,  indicating  that  it  may  be  applicable  to 
arbitrary  wave  fields,  both  linear  and  nonlinear. 
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APPENDIX  MATLAB  CODE 


This  appendix  contains  the  particular  mfile  (Matlab  Code)  that  we  used  to 
analyze  the  data  from  the  three  buoys  in  order  to  extract  the  information  that  we 
were  interesting  on. 

%Loading  the  desired  data 
load  16jul1a.dat 
p1=16jul1a(:,1) 
clear  16jul1a 
load  16jul2a.dat 
p2=jul2a(:,1) 
clear  16jul2a 
load  16jul3a.dat 
p3=16jul3a(:,3) 
clear  16jul3a 
p1=p1-mean(p1) 
p2=p2-mean(p2) 
p3=p3-mean(p3) 
fs=2 

M1=  max(size(p1)) 

N=1 72800 
f=0:fs/N:(N-1)*(fs/N) 
tbp=0.5 

t1=[1:M1].*(tbp/60) 

%  Looking  for  the  FFT  of  the  data  from  the  three  buoys  to  select  the 
desired  frequencies 

P1=fft(p1)*2/length(p1)  ; 

P2=fft(p2)*2/length(p2)  ; 

P3=fft(p3)*2/length(p3)  ; 

%  Applying  the  filter  and  Hilbert  Transform  to  the  data 

cfjow=0.08  ;  %cut  off  frequency  for  buoys  1 ,2  and  3. 

[b1,a1]=cheby2(5,40,cfJow/(fs/2));%Chebyshev  filter  5  order  40  dB  for  data 

%form  buoy  1 


%Sampling  frequency 
%Get  Number  of  data  points 
%Total  number  of  points 
%Normalizing  the  x-axis  to  frequency 
%Time  between  data  points  in  seconds 
%Normalizing  the  x-axis  to  time  in  minutes 
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Y1  =filtfilt(b1 ,  a1 ,  p1 )  ;  %Filtering  data  from  buoy  1 

clear  p1  ;  %Clear  data  pressure  1  from  memory 

V_wavep1=hilbert(Y1)  ;  %Hilbert  Transform  data  from  buoy  1 

phasep1=angle(V_wavep1)  ;  %Phase  of  the  Hilbert  Transform  for  buoy  1 
[S1  ,f1]=psd(Y1 ,1024,fs)  ;%examine  spectrum  of  the  low  wave 

%without  filtering 

[S2,f2]=psd(V_wavep1 ,1 024,fs)  ;%examine  spectrum  of  the  low  wave 

%after  filtering  and  the  Hilbert  transform 
[b2,a2]=cheby2(5,40,cf_low/(fs/2));  %Chebyshev  filter  order  5 , 40  dB. 
Y2=filtfilt(b2,  a2,  p2)  ;%Filtering  data  from  buoy2 

clear  p2  ;  %Clear  data  pressure  2  from  memory 

V_wavep2=hilbert(Y2)  ;  %Hilbert  Transform  data  from  buoy  2 

phasep2=angle(V_wavep2)  ;  %Phase  of  the  Hilbert  Transform  for  buoy  1 
[b3,a3]=cheby2(5,40,cfJow/(fs/2));  %Chebyshev  filter  for  data  from  buoy  3 
Y3=filtfilt(b3,  a3,  p3)  ;  %Filtering  data  from  buoy  3 

V_wavep3=hilbert(Y3)  ;  %Hilbei1  Transform  data  from  buoy  3 

phasep3=angle(V_wavep3)  ;  %Phase  of  the  Hilbert  Transform  for  buoy  t 
clear  p3  ;  %Clear  data  pressure  3  from  memory 

%Looking  in  a  two  hours  windows  of  the  data  for  each  buoy 
V_wavep1  =V_wavep1  (1:1 4400)  ;  %points  correspond  to  a  two  hours  windows 
V_wavep2=V_wavep2(1 : 14400)  ;  %points  correspond  to  a  two  hours  windows 
V_wavep3=V_wave  p3(1 : 14400)  ;  %points  correspond  to  a  two  hours  windows 
Y1=Y1  (1 : 14400)  ;  %points  correspond  to  a  two  hours  windows 

Y2=Y2(1 : 14400)  ;  %points  correspond  to  a  two  hours  windows 

Y3=Y3(1 ;  14400)  ;  %points  correspond  to  a  two  hours  windows 

t1=t1  (1 :14400)  ;  %Two  hour  time  window 

%Looking  for  the  phase  of  the  Hilbert  Transform  for  the  two  hour  window 
theta1=unwrap(angle(V_wavep1));  %Phase  of  the  Hilbert  transform  in  radians 
theta2=unwrap(angle(V_wavep2));  %Phase  of  the  Hilbert  transform  in  radians 
theta3=unwrap(angle(V_wavep3));  %Phase  of  the  Hilbert  transform  in  radians 
%Looking  for  the  derivative  of  the  phase  of  the  Hilbert  Transform  for  the 
two  hour  window 

dthetal  =  diff(theta1  ./tbp)  ;  %By  definition  the  Local  frequency  at 
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%buoy  1 

dtheta2  =  diff(theta2./tbp)  ;%By  definition  the  Local  frequency  at 

%buoy  2 

dthetaS  =  diff(theta3./tbp)  ;%By  definition  the  Local  frequency  at 

%buoy  3 

%Looking  for  the  cross-correlation  of  the  phase  of  the  Hilbert  transform 
Cp21=xcorr(phasep2,  phasepl.'coeff)  ;%Xcorrelations  of  the  phase  of  the 

%signal  from  Buoy  2  and  1 . 

Cp31=xcorr(phasep3,  phasep1,'coeff')  ;%Xcorrelations  of  the  phase  of  the 

%signal  from  Buoy  3  and  1 . 

Cp32=xcorr(phasep3,  phasep2,'coeff')  ;%Xcorrelations  of  the  phase  of  the 

%signal  from  Buoy  3  and  2. 

Cps21  =fftshift(Cp21 )  ;  %Shift  of  the  correlations 

Cps31  =fftshift(Cp31 )  ;  %Shift  of  the  correlations 

Cps32=fftshift(Cp32)  ;  %Shift  of  the  correlations 

%Plots  of  the  FFT  of  the  data  pressure  from  the  three  buoys 
subplot  (3,1,1) 

plot(f,abs(P1)),axis([0  0.5  0  0.05]),grid,title('FFT  of  Data  Pressure 
of  buoy  1  vs  Frequency'), 
subplot  (3,1,2) 

plot(f,abs(P2)),axis([0  0.5  0  0.05]), grid, 
title('FFT  of  Data  Pressure  of  buoy  2  vs  Frequency'), 


33 


subplot  (3,1,3) 

plot(f,abs(P3)),axis([0  0.5  0  0.05]),grid, 

title('FFT  of  Data  Pressure  of  buoy  3  vs  Frequency'), 

xlabel('Frequency'),ylabel('Normalized  magnitude')  ; 

clear  PI ;  clear  P2;  clear  P3  ;%Release  memory 

%Plot  of  the  Power  Estimate  Spectrum 

figure;  clf; 

subplot(2, 1,1) 

loglog(f1,S1),axis([1E-3  1  IE-6  1]); 

title('Power  spectrum  Estimate  of  Pressure  data  frequency  without  filtering'); 
ylabel('magnitude  (dB)'); 
subplot(2,  1 , 2) 

loglog(f2,S2),axis([1E-3  1  IE-6  1]); 

title('Power  spectrum  Estimate  of  V_wavep1  after  Hilbert  Transform  '); 
xlabel(['frequency(Hz)-cutoff  frequency',num2str(cfJow),'Hz']); 
ylabel('magnitude  (dB)'); 

%Plot  of  the  Phase  in  radians  (Unwrap)  of  the  Hilbert  transform 

figure;  clf; 
subplot(3,  1,1) 
plot(t1,(phase_1)), 

title('Plot  of  the  phase  of  Hilbert  Transform  of  the  pressure  data  buoy  T); 
xlabel('time(min)');  ylabel('Phase(rad)'); 
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subplot(3,  1,2) 

plot(t1,(phase_2)),title('Plot  of  the  phase  of  Hilbert  Transform  of  the  pressure 
data  buoy  2');  xlabel('time(min)');ylabel('Phase(rad)'); 
subplot(3,  1 , 3) 

plot(t1,(phase_3)),title('Plot  of  the  phase  of  Hilbert  Transform  of  the  pressure 
data  buoy  3');  xlabel('time(min)');  ylabel('Phase(rad)'); 

%Plots  of  the  Local  frequency  as  a  function  of  the  time,  data  from  buoy  1 

figure;  clf; 

plot(t1,dtheta1);  title('Plot  of  the  Local  Frequency  as  a  function  of  time  at  buoy 
T);  xlabel('time(min)');  ylabel('frequency(rad/sec)'); 

%Plots  of  the  correlations  of  the  data  of  the  three  buoys 

figure;clf; 
subplot(3,  1,1) 

plot(t1,abs(Cp21)),grid,title('Correlation  of  Phase  of  the  data  from  Buoy  2  and 
Buoy  1'),xlabel(Time(sec)') ; 
subplot(3,  1 , 2) 

plot(t1,abs(Cp31)),grid,title('Correlation  of  Phase  of  the  data  from  Buoy  3  and 
Buoy  1'),xlabel(Time(sec)'); 
subplot(3,  1 , 3) 

plot(t1,abs(Cp32)),grid,title('Correlation  of  Phase  of  the  data  from  Buoy  3  and 
Buoy  2'),xlabel(Time(sec)'); 

%Pi6t  of  the  shift  of  the  correlations  of  the  Phase  from  the  three  buoys 
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figure;clf; 
subplot(3,  1,1) 

plot(t1  ,abs(Cps21)),grid,title('Shift  of  Correlation  of  Phase  of  the  data  from  Buoy 

2  and  Buoy  1'),axis([0  1500  0  1]); 
subplot(3, 1 , 2) 

plot(t1,abs(Cps31)),grid,title('Shift  of  Correlation  of  Phase  of  the  data  from  Buoy 

3  and  Buoy  1'),axis([0  1500  0  1]); 
subplot(3,  1 , 3) 

plot(t1,abs(Cps32)),grid,title('Shift  of  Correlation  of  Phase  of  the  data  from  Buoy 
3  and  Buoy  2'),axis([0  1500  0  1]); 

%Plot  of  the  Maximum  values  of  the  shifted  correlations 

figure:clf; 
subplot{3,  1,1) 

plot(abs(Cps21)),grid,title('Max.  value  of  the  Correlation  of  the  Phase  of  the  data 
from  Buoy  2  and  Buoy  1'),axis([0  40  0.7  0.85]); 
subplot(3,  1 , 2) 

plot(abs(Cps31)),grid,title('Max.  value  of  the  Correlation  of  the  Phase  of  the  data 
from  Buoy  3  and  Buoy  1'),axis([75  150  0.5  0.7]); 
subplot(3,  1 , 3) 

plot(abs(Cps32)),grid,title('Max.  value  of  the  Correlation  of  the  Phase  of  the  data 
from  Buoy  3  and  Buoy  2'),axis([120  200  0.5  0.66]); 
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